
function [c,ceq,Gc,Gceq]=eq_con_jac_NOlearn(pars,~,~,M,X,Z,sgi_nodes,sgi_weights,ep_sd,~)

c=[];

[N,T,K]=size(X);
KZ=size(Z,3);

M=reshape(M,N*T,1);

alpha=repmat(pars(1:N,1),T,1);
theta=pars(N+(1:2),1);
beta=pars(N+2+(1:2*N),1);
v=pars(3*N+2+(1:N),1);
f=repmat(pars(4*N+2+(1:N),1),T,1);
% vs=pars(5*N+2+(1:N),1);
xi=pars(6*N+2+(1:K),1);
gamma=pars(6*N+2+K+(1:KZ),1);
kappa=pars((6*N+2+K+KZ+1):end,1);

ep_sd=repmat(ep_sd,T,1);

if nargout<=2
    
    ZZ=zeros(N,N);
    for kz=1:KZ
        ZZ=ZZ+Z(:,:,kz)*gamma(kz,1);
    end 
    Q=diag(v.*ep_sd(1:N,1))*exp(-ZZ)*diag(v.*ep_sd(1:N,1));
    P=kron(eye(2),tril(Q)+tril(Q,-1)');
    P=P\eye(2*N);
    
    XX=zeros(N,T);
    Pd=repmat(sqrt(diag(P\eye(2*N))),1,T);
    B=repmat(beta,1,T);
    for t=1:T
        XX(:,t)=permute(X(:,t,:),[1,3,2])*xi;
    end
    XX=reshape(XX,N*T,1);
    PdA=reshape(Pd(1:N,:),N*T,1);
    Pd=reshape(Pd(N+(1:N),:),N*T,1);
    BA=reshape(B(1:N,:),N*T,1);
    B=reshape(B(N+(1:N),:),N*T,1);

    UU0=exp(alpha+theta(1)*(PdA*sgi_nodes(:,1)'+BA+ep_sd*sgi_nodes(:,2)'));
    UU0=M.*(UU0./(1+UU0));
    UU1=exp(alpha+theta(2)*(Pd*sgi_nodes(:,1)'+B+ep_sd*sgi_nodes(:,2)')-f-XX-kappa);
    UU1=M.*(UU1./(1+UU1));

    ceq=(UU1-UU0)*sgi_weights;
    % ceq=[U1,U0];
    
else
    
    ZZ=zeros(N,N);
    for kz=1:KZ
        ZZ=ZZ+Z(:,:,kz)*gamma(kz);
    end 
    Q=diag(v.*ep_sd(1:N,1))*exp(-ZZ)*diag(v.*ep_sd(1:N,1));
    P=kron(eye(2),tril(Q)+tril(Q,-1)');
    P=P\eye(2*N);
    P=repmat(P,1,1,T);

    XX=zeros(N,T);
    Pd=repmat(sqrt(diag(P(:,:,1)\eye(2*N))),1,T);
    B=repmat(beta,1,T);
    for t=1:T
        XX(:,t)=permute(X(:,t,:),[1,3,2])*xi;
    end
    XX=reshape(XX,N*T,1);
    PdA=reshape(Pd(1:N,:),N*T,1);
    Pd=reshape(Pd(N+(1:N),:),N*T,1);
    BA=reshape(B(1:N,:),N*T,1);
    BD=reshape(B(N+(1:N),:),N*T,1);

    yA=PdA*sgi_nodes(:,1)'+BA+ep_sd*sgi_nodes(:,2)';
    yD=Pd*sgi_nodes(:,1)'+BD+ep_sd*sgi_nodes(:,2)';
    
    UU0=exp(alpha+theta(1)*yA);
    UU0=M.*(UU0./(1+UU0));
    UU1=exp(alpha+theta(2)*yD-f-XX-kappa);
    UU1=M.*(UU1./(1+UU1));

    ceq=(UU1-UU0)*sgi_weights;

    UU0=UU0.*(1-UU0);
    UU1=UU1.*(1-UU1);

    Gc=[];
    
    Gceq=zeros(N*T,length(pars));
    % alpha
    Gceq(:,1:N)=((UU1-UU0)*sgi_weights).*repmat(eye(N),T,1);
    % theta
    Gceq(:,N+1)=-(UU0.*yA)*sgi_weights;
    Gceq(:,N+2)=(UU1.*yD)*sgi_weights;
    % betaA
    Gceq(:,N+2+(1:N))=(-theta(1)*UU0*sgi_weights).*repmat(eye(N),T,1);
    % betaD
    Gceq(:,N+2+N+(1:N))=(theta(2)*UU1*sgi_weights).*repmat(eye(N),T,1);
    % v, gamma
    DS=zeros(2*N,N+1);
    for n=1:N
        DS(:,n+1)=DB_vg_NOlearn(Z,ZZ,v,ep_sd(1:N,1),n);
          % v
        Gceq(:,3*N+2+n)=(theta(2)*UU1.*(repmat(DS(N+(1:N),n+1),T,1)*sgi_nodes(:,1)')-...
        theta(1)*UU0.*(repmat(DS(1:N,n+1),T,1)*sgi_nodes(:,1)'))*sgi_weights;
    end
    DS(:,1)=DB_vg_NOlearn(Z,ZZ,v,ep_sd(1:N,1),0);
    Gceq(:,6*N+2+K+(1:KZ))=(theta(2)*UU1.*(repmat(DS(N+(1:N),1),T,1)*sgi_nodes(:,1)')-...
        theta(1)*UU0.*(repmat(DS(1:N,1),T,1)*sgi_nodes(:,1)'))*sgi_weights;
    % f
    A=-(UU1)*sgi_weights;
    Gceq(:,4*N+2+(1:N))=A.*repmat(eye(N),T,1);
    % kappa
    Gceq(:,(6*N+2+K+KZ+1):end)=diag(A);
    % xi
    for k=1:K
        Gceq(:,6*N+2+k)=-(UU1.*reshape(X(:,:,k),N*T,1))*sgi_weights;
    end
    
    Gceq=sparse(Gceq');
    
end


